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ABSTRACT 

We consider the influence of potential quark matter existing at high densities in neu- 
tron star interiors on gravitational waves (GW) emitted in a binary neutron star 
merger event. Two types of equations of state (EoS) at zero temperatures are used, 
one describing pure nuclear matter, the other nuclear matter with a phase transition 
to quark matter at very high densities. Binary equilibrium sequences close to the in- 
nermost stable circular orbit (ISCO) are calculated to determine the GW frequencies 
just before merger. It is found that EoS effects begin to play a role for gravitational 
masses larger than M^o — 1-5Mq. The difference in the gravitational wave frequency 
at the ISCO grows to up to ~ 10% for the maximal allowed mass given by the EoSs 
used. Then, we perform 3D hydrodynamic simulations for each EoS varying the initial 
mass and determine the characteristic GW frequencies of the merger remnants. The 
implications of quark matter show up mainly in a different collapse behaviour of the 
merger remnant. If the collapse does not take place immediately after merger, we find 
a phase difference between two EoS's in the post-merger GW signal. We also compare 
the GW frequencies emitted by the merger remnant to values from simulations using 
a polytropic EoS and find an imprint of the non-constant adiabatic index of our EoSs. 
All calculations are based on the conformally flat (CF) approximation to general rel- 
ativity and the GW signal from the merger simulation is extracted up to quadrupole 
order. 

Key words: hydrodynamics, gravitational waves, equation of state, stars: neutron 



1 INTRODUCTION 

Binary neutron star mergers belong to the strongest 
gravitational wave sources for interferometer-type gravita- 
tional w ave observatories as LIG O (iAbramo yici et al.,.19921 . 
VIRGO fBradaschia et al.lll990D . GEO600 jLucklll997H and 



TAMA jAndo et al.l l200ll) . After a long inspiral process 
which lasts millions of years, the final merger phase takes 
place on a millisecond timescale. The onset of merger is 
given when the two companions become dynamically unsta- 
ble near the innermost stable circular orbit (ISCO) and mass 
transfer starts. Since the process has no definite symmetries, 
3D hydrodynamic simulations are necessary. In addition, be- 
cause of the compactness of the system, general relativistic 
effects have to be considered and treated as accurately as 
possible. 

It has been pointed out that the gravitational wave fre- 
quency few relates to the compactness {M/R)oo of a neu- 
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tron star as /gw ~ (M/R)^^ near the ISCO, meaning 
that the gravitational wave frequency just before merger 
carries information on the radius of the neutron star. 
Hence the GW signal may constrain the e quation of state 
(EoS) of high density matter, (see e.g. iLai fc WisemanI 
|l996; Rasio & Shapiro, J_999) • It has also been found re- 
cently that a neutron star formed after the merger may 
be supported by rapid and differential rotation even if its 
mass exceeds 60% of the maximum mass of a single non- 
rotating neutron star, and the GW are emitted due to non- 
axisymmetric and quasiradial oscillat ions of the remnant for 
longer than the dynamical timescale Ishibata fc Urvull2000l . 
120021: IShibata. Taniguchi fc Urvul2003tl . If these oscillations 
persist for a longer time, the integrated gravitational wave 
may be detectable, carrying information on the high density 
matter. 

In addition to general relativistic gravity, various 
physics should be taken into account for many aspects of 
the binary neutron star merger problem such are the nuclear 
forces summarized in the EoS, neutrino physics and mag- 
netic fields. So far, however, due to this complexity, investi- 
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gations have concentrated either on the relativistic aspects 
of the problem using a simple EoS, or on the microphys- 
ical aspects while treating gravity in a Newtonian frame- 
wor k. The pioneering work was d one by Oohara & Nakamura 
in llOohara fc Nakamural Il999l and references therein). 
Relativistic aspects have been considered using the post- 
Newt onian approximation jAval et alJl200ll : iFaber fc Rasid 
|2002| . and references therein), the conformally flat ap- 
proximation ( Mathews fc Wilsoru.2000. : , Oechslin et al. 2002; 
iFaber. Grandc lcmcnt fc Rasidl2003l and re ferences therein) 
and a fully general relativi stic treatment Jshibata fc Urvul 
I2OOOI. I2OO2 : Shibata, Tani euchi fc Urviil 12001 and refer- 
ences t herein). Microphysical i mprovements have been ap- 
plied by|RAiffert_fc_Jaiika 1I2OOIL and references therein) and 
iRosswog fc Dav iea ll200 2. and references therein) using dif- 
ferent EoSs llShen et alJll998tlLattimer fc Swestvll99lfl and 
a leakage scheme to account for neutrino e mission after the 
^mer^e r. For a review on the topic, see llRasio fc Shapirol 

In this paper, we focus on two EoS for high density 
neutron star matter, one describing pure nuclear matter, 
the other nuclear matter with a transition to a quark phase 
at very high densities. A reason for investigating these EoSs 
is that the quark phase transition may be one of the most 
dramatic phenomena to change the compactness of neutron 
stars, and therefore its influence may be clearly observed in 
the GW spectrum of an inspiralling binary system. Following 
our previous work dOechslin et al.ll2002h . we consider the 
merger problem in the conformally flat approximation. 

The paper is organised as follows. In sect.|21 we summa- 
rize the formalism to solve the Einstein and the relativistic 
hydrodynamic equations, and describe the numerical imple- 
mentation, the choice of EoS and the initial conditions. In 
sect. 13 we present the results of the paper and flnally, we 
draw conclusions in sect. |1| 



2 FUNDAMENTAL INGREDIENTS 

In this section we describe the numerical methods and the 
physics on which our simulations are based. We consider a 
general relativistic fluid whose internal properties are de- 
scribed by a given EoS. Two sets of equations have to be 
solved simultaneously. On one hand side, we consider the 
relativistic hydrodynamic equations governing the fluid mo- 
tion on the other side the Einstein equation of general rel- 
ativity determining the spacetime metric and therefore the 
gravitational interaction. On top of it, an EoS which closes 
the system of hydrodynamic equations, has to given as an 
input. 

2.1 Relativistic Hydrodynamics and 

Approximation for General Relativistic 
Gravity 

In the 3+1 decomposition of spacetime, the metric ds^ — 
gfj,vdx^dx^ can be written as 

da^ = {-a + /3,f3')dt'^ + 2IJidx'dt + -y^jdx'dx^ (1) 

where a is the lapse function, /3' is the shift vector, and 
is the spatial metric. Spacetime quantities are decomposed 
with respect to the foliation using the hypersurface normal 



deflned as = —ad^t with n'^n^ = — 1 and a projection 
tensor defined as 7^,^ = g^v+n^rii,, whose spatial component 
agrees with the spatial metric. 

The Einstein field equations can be written as a set of 
two evolution equations 

^t'y^J = ~2aK,j + + "^^(3,, (2) 
~9,^lS^J+A^l-^ij(S ~ pe)] 

(3) 

and two constraint equations 

R + K^-K,,K'' = WnpE, (4) 
\7j{K'' -j'^K) = 8Tvf (5) 

for the dynamic variables 7^-, and Kg, the extri nsic curva- 
ture of the hypersurface (see e.g. iBaumgarte a nd ShaDi3 
|200Sj). Here, Rij is the Ricci tensor and V the covariant 
derivative associated with 7^^ . The stress energy tensor T^i^ 
is decomposed into pE = n^n'^T^,^, the matter energy den- 
sity, = 7^n^r^'', the matter momentum density and 
Sij — 71^7^1/7^", the spatial projection of the stress energy 
tensor. Finally, R, S and K denote the traces of Rij, Sij and 
Kij, respectively. 

In the following discussion, we consider a perfect fiuid 
with a matter stress energy tensor 

Tfj,„ = phufj,u„ + pg^i, (6) 

Here, p refers to the rest mass density, h = l-\-p/ p + e to the 
specific relativistic enthalpy, e to the specific internal energy, 
to the four velocity and p to the fluid pressure. Then, 

pE = ph{au')^-p (7) 
f = phau°u^'Yl. (8) 

The Lorentz factor W — au^ can be calculated using the 
normalisation condition Upit^ = —1. 

au'' = {l + j'u,Uj)'-/\ (9) 

Isenberg proposed a waveless approximation to gen- 
eral relativity, in which he truncates some terms in Ein- 
stein's equation written in th e ADM forma lism to de- 
duce an elliptic type formalism lllsenbe rg"l978') . Later, the 
same set of eq uations have been rediscovered by W ilson 
and Mathews dWilson. Mathews fc Marronettil Il99d) and 
widely us ed to solve si ngle or binary neutron stars prob- 
lems. (e.gJOechslh^t^L 2002) as well as binary black hol e 
systems l|Grandclement. Gourgoulhon fc Bonazzolal l2002ll . 
The Isenberg- Wilson-Mathews theory has its own Hamil- 
tonian as pointed out by Isenber g himself (also discussed in 
iFriedman. Urvti fc Shibatall2002ll . 

In this approximation, the conformally flat condition 
for the spatial geometry "fij — ip'^Sij and K = dtK = are 
imposed, where 5ij is the flat 3-metric. It may be viewed 
that the second one is a choice of a temporal gauge condi- 
tion, the maximal slicing condition, while the first one is an 
approximation as w ell as partly a spatial g auge choice (de- 
tails can be found in lBaumgarte et alJll99q) . The remaining 
metric variables (V',Q,/3*) do not satisfy all components of 
the Einstein's field equation consistently. We pick the con- 
straints and the trace of the evolution equation for Kij as 
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equations for the above five variables, which leads to the 
same set of equations derived from Isenberg's Ifamiltonian 
mentioned above. 

The approximation leads to a considerable simplifica- 
tion of the Einstein equations, since all metric equations 
can be written in elliptic form. 

The trace of the evolution equation for Kij together 
with the maximal slicing condition K — leads to an equa- 
tion for the lapse 

A(aV) = 2TTa^^{pE + 2S) + laiP^K^jK'^, (10) 

o 

The trace-free part of the evolution equation for 
together with the conformally flat condition 7ij — 'tp'^Sij 
provides 

2ail,-^K,j = djPr + drk ~ |<5^J9fc/3^ (11) 

where di is a derivative with respect to the coordinate asso- 
ciated with 5ij and Pi relate to the shift through the flat 
metric as Pi = SijP-'. The Hamiltonian constraint provides 
an equation for the conformal factor 



where 



(12) 



Finally, with the momentum constraint we obtain an equa- 
tion for the shift vector 



A/3' + -d'djP' = dj In 



+ 167rQ7/)*j', 



where 9' = S'^dj. Using the definition 



P'=B^- ^d'x, 
I13II splits into two simpler parts 



(13) 



(14) 



and 



(15) 



Ax = d,B' 



(16) 



While the trace-free components of the evolution-equation 
for @ are used to relate Kij and the metric variables 
as in equation lllll . its trace part and the trace- free com- 
ponents of the evolution-equation for Kij Q are dropped 
by the CF approximation. The trace of ^ can be used to 
check the accuracy and reliability of the CF approximation 
JPimmelmeier et aLlliool lOechslinlliooS) 

The relativistic hydrodynamic equations can be written 
in the following way 



dtp* + di(p*v') 
dt{p*e) + 9i(p*eu') 



0, 



a I 1/2n 
dt\OLU 7 ' ) 



dt{p*huk) + di{p*hukv^ 
p*h 



I a / 1/2 i\ 

-^-Oi[au 7 ' ti ) 



a 



2uO 



■9fc7*' 



(17) 

(18) 
(19) 
(20) 



pau det(7ij) and 



-P' + 



In the case of conformal fiatness, the system reduces to 
d , 

^aip^dip — aifdiO + UjdiP^ 



-p'div' 



d 



p' 

2ufcUfc 



dt 
with 

Ui = hui 



* 0,6 
p = pau tp 



— OiV Otintau Tp ), 

p* p* 



and 



(21) 
(22) 

(23) 

(24) 
(25) 



(26) 
(27) 

(28) 



Here, we have cast the system into a Lagrangian formula- 
tion which is better suited for an implementation on a com- 
puter using a Lagrangian scheme like SPH (cf. sect. l2.2(l . The 
conserved hydrodynamic variables are {p* , Ui, e) whereas the 
primitive variables are (p, v',e). 



2.2 Numerical implementation 

To solve the system of hydrodynamic equations 1231 - 
1251 1 we use the Smoothed Particle Hydrod ynamics method 
(SPH; jBendll990l: iMonaehan fc Gingoldl ll983')'). which is 
widely used in astronomical and astrophysical simulations. 
The original form of SPH that solves the Newto nian hy- 
drodynamic equations is presented in llBenj|l990l) . We list 
here the generalized version which is appropriate for the 
above relativistic hyd rodynamic equations. For details, see 
('Oe chsUn et alj|2002ri . The continuity equation turns into a 
relation between p* and the rest masses of the individual 
particles 



Pa = ^ mtWab, 



(29) 



where mj, is the rest mass of particle b and Wab = W^(|ra — 
ffcl, hspn) denotes the weight given by the standard spheri- 
cal spline Kernel function W{r,hspH)- The rest masses mt 
and the smoothing length hspH are initially chosen to fit 
an initial spatial distribution of p*. The pressure gradient 
contained in the momentum equation 12511 is calculated as 
in standard SPH but replacing p by p*, 

\^^pa = ~ymb + ^] d.Wat. (30) 

The energy equation 12511 contains a pdV-term similar to 
Newtonian SPH 



P 



di' 



IE 



nib 



Pa ^ Pb 
Pa Pa Plpb 



{vl - vl)d,Wab. 



(31) 

The additional terms in the momentum and energy equation 
that arise from the gravitational interaction are evaluated by 
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first solving the equations for the metric variables a,-;/) and 
/3' on a overlaid grid (see discussion below) and by mapping 
back those quantities onto the particles by second-order in- 
terpolation. The total time derivative of \n{avP^^)a in the 
energy equation is evaluated using second order finite differ- 
encing in time. 

On top of these equations, an artificial viscosity (AV) 
scheme is implemented wit h time-dependent viscosity pa- 
rameters iMorris fc Monaghan 1997) which have only sig- 
nificant values in the presence of shocks. The AV produces 
an additional viscous pressure which is added to the physical 
fluid pressure (Sioglcr & Riffcrt 2000). 

The field equations 110ll2ll5l and [T^ are discretized on 
a grid which covers the matter distribution. The evaluation 
proceeds generally in three steps 

• First, we calculate the source corresponding to the po- 
tential. The involved hydrodynamic quantities defined on 
the SPH particles are transferred onto the grid by assigning 
SPH-interpolated values to the grid points, i.e. 

{f{x))sPH = /^■^W'd^ " r,!, ft,). (32) 

• Afterwards, the potential is obtained by solving the cor- 
responding Poisson-equation. The ip- and the at/i-equation 
are closely linked together and are therefore solved in two it- 
eration steps using the solution from the previous timestep 
as a guess. The remaining four equations determining the 
shift vector are also closely coupled and are solved simulta- 
neously within the FMG algorithm. 

• All needed derivatives are calculated on the grid using 
finite differencing. 

• Finally, the potential and its derivatives are mapped 
back onto the SPH particle dis tribution using a triangu - 
lar shaped cloud (TSC) method jHockev fc EastwooJl992l) 
known from particle mesh codes. This is equivalent to a sec- 
ond order interpolation. 

Boundary conditions and extensions of the solution be- 
yond the grid are obtained using a multi pole expansion 
of th e source term. For details we refer to iOechslin et alJ 
|2003). 

Like the Newtonian and the first-order post-Newtonian 
(IPN) approximation, the CF approximation does not in- 
clude gravitational radiation and its reaction by construc- 
tion. We therefore have to add an additional gravitational 
wave (GW) extraction scheme and a radiation reaction force 
which accounts for the angular momentum and energy car- 
ried away by GW. The waveform in trans verse-traceless 
gauge is extracted in the slow-m otion limit jThornelllQSOl: 
IWilson. Mat hews fc Marronettll99&'l and up to quadrupole 
order using 

hf/'{t,r) = ^PijkiQki{t-r), (33) 

where 

Q^, = STF |-2 J S^{x)xix.id^x^ (34) 

is the mass-quadrupole of the system. The radiation-reaction 
force is chosen in a similar way to the original Burke- Thorne 
expression = adiV with V = —l/5xiXjQl'^\ but its 
strength a is chosen such that the resulting angular mo- 

. • tot • num. 

mentum loss by backreaction Ji — Ji reproduces the 



expression 

j. = -|6,,feQfiQi!i (35) 

of the slow motion GW extraction approximation up to a 
difference of the order of ~ 5%. Since Ji = eij^ J p*{x^Uk — 
x^Uj) (f'x, we can obtain the angular momentum change us- 
ing 

Ji = eijk j p* (^{x^Uk — x'^iij) + {v^Uk - i]''uj)^ d^x. (36) 

The expressions j/^^"^ , the numerical angular momentum 
error, and Ji ° , the total angular momentum change per 
unit time, are obtained by evaluating (1361 before adding 
and after adding the backreaction force terms, respectively. 
Since Ji has an oscillating behaviour around zero, as 
a result thereof, Ji ° depends on the binary orbit. There- 
fore, in practical determination of the value of cr, we take 
an average {Ji ° — Ji""^) for about a half orbit, and com- 
pare with (J) computed from the quadrupole formula 1351 
that is averaged for the same orbit. For model Bl we obtain 
for the ratio (Ji "'"")/( ~ 0.3% at the ISCO, where 
the average is made over one orbital period, while at a cer- 

. • num • tot 

tarn period of time Ji may become comparable to Ji 
Since the ratio [Ji ) / \Ji ) is negligible, the radiation re- 
action force is clearly driving the inspiral process, although 
the predicted time to merger have to be taken with care 
due to the approximative radiation reaction scheme and the 
spatially conformal flat assumption used in our initial condi- 
tions and simulations. Fully rela tivistic inspiral simulations 
iMiller. Gressman fc SuenI (|200^ are intrinsically consistent 
concerning radiation reaction dissipation but they still de- 
pend on the outer boundary location. 

2.3 Equation of state 

To investigate the influence of quark matter to a neutron 
star merger, we consider two equations of state, an EoS de- 
scribing pure nuclear matter ("hadronic EoS") and an EoS 
describing nuclear matter with a phase transition to quark 
matter at very high densities ( "hybrid EoS" ) . 

The hadronic EoS is realized using the non-linear a-uj- 
model in the relativistic mean field approximation with the 
TMl parameter set which is motivated by a least-square fit 
to experimental results including stable and unstable nuclei 
(Suga hara fc Toki 1994). At densities above p > 10^''g/cm^ 
this is a good approximation. For lower values, when inho- 
mogeneous nuclear matter appears and the non-linear a-uj 
model is not valid any more, ot her approximations have to 
be used (e.g. IShen et al.l 1199^1 . In our case, we append a 
polytropic EoS p = Kff where k and F are adjusted to fulfill 
smoothness of pressure and internal energy. We choose the 
transition density to the polytropic EoS at 2 x 10^'*g/cm'^ 
which leads to P ~ 2.86, similar to common realistic EoSs in 
this density regime. The hybrid EoS is obtain ed by combin- 
ing the hadronic EoS with a MIT bag model iChodos et alJ 
1 974 

using a v ariable pressure phase transition construction 
liGleiidenninslll992ft . We use massless up and down quarks 
and a Bag constant of 90MeV/fm'^. The resulting hybrid 
EoS then describes three physical phases 

• A pure hadronic phase below ~ 5 x lO^^'g/cm^ ~ 1.8po 
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Figure 1. Comparison of the considered nuclear EoS. The hybrid 
EoS (dashed curve) has a phase transition region (1.8po < P < 
3.5po) where adiabatic indices are substantially lower, followed 
by a quark phase which is has a similar stiffness as the nuclear 
matter at those densities. 



where the EoS coincides with the hadronic EoS. Here, po := 
2.8 X 10^'*g/cm^ is the nuclear saturation density. The stiff- 
ness in this region varies between F « 3 and F ~ 2.5. 
• A quark-hadron mixed phase between 



5 X 10"g/cm^ 



and ~ 10 g/cm ~ 3.5po where both quark and hadrons 
are present. In this phase transition region, the EoS sub- 
stantially softens with F « 1 — 1.5. 

• A pure quark phase above « lO^^g/cm^ where the MIT- 
bag-model EoS with similar adiabatic index but lower abso- 
lute pressure in comparison to the hadronic EoS is applied. 

In both EoSs, we neglect temperature effects, i.e. we set 
T=0. The redundant internal energy information from the 
EoS is dropped and we consider the pressure as a function 
of the density alone. This is a good approximation in the 
high-density regime where thermal effects to the pressure 
are rather small. At lower densities (e.g. in the disk around 
the merger remnant) the thermal component is certainly not 
negligible. However we concentrate in this work on the high- 
density merger remnant so we claim to be rather accurate 
with this approximation although thermal effects could have 
an effect via neutrino losses. 



2.4 Initial Conditions 

To construct initial conditions, we assume the binary to be 
in a quasi-equilibrium state. This assumption is not very 
well satisfied in the very vicinity of the ISCO ( Miller 2003) 
and we therefore choose slightly wider initial orbits. In our 
models the binaries have initially no radial velocity either 
and need about half an orbit to begin with the inspiral. We 
apply the method of Uryri & Eriguchi (2000) to construct an 
initial configuration. This method takes an EoS as an input 
and then solves the hydrostatic equation for an irrotational 
velocity field together with the Einstein field equations in 
the CP approximation. A similar method has b een devel- 
oped by (iBonazzola. Gourgoulhon fc Marcklll999^ . We map 



Table 1. Input parameters of the considered models using the 
hadronic and hybrid EoS. Mao denotes the gravitational mass of 
one NS in isolation and C is the compactness. Both the masses 
and the compactness refer to one binary component in isolation. 
Nparticle denotes the number of SPH particles and NQ^id is the 
number of grid points. In model El, 97^ gridpoints are used, 
but at the onset of collapse, the resolution is enhanced to 129^. 
Models with a '1' in the label are calculated with the hadronic 
EoS, models with a '2' used the hybrid EoS. 



Model 


MooIMq] 


C = (M/R)ac 


^ particle 




Al 


1.2 


0.1276 


154038 


973 


B1,B2 


1.35 


0.1424 


153993 


973 


C1,C2 


1.4 


0.1475 


169798 


973 


Dl 


1.5 


0.1577 


100488 


973 


D2 


1.5 


0.1579 


100548 


973 


El 


1.75 


0.1853 


162536 


973(1293) 


E2 


1.75 


0.1959 


157222 


1293 



Table 2. Selected properties of our initial models resulting from 
the input parameters shown in TablelTl Mq is the restmass of one 
single NS, P the orbital period, fcwfi the corresponding GW 
frequency, do the initial orbital separation and Roc, the radius 
of one single isolated NS measured in Schwarzschild coordinates. 
Note that the GW frequencies are lower than those in Fig.l2lwhich 
are taken at the ISCO. 



Model 


Mo[Mq] 


P[ms\ 


faw.o = 2/P 


do [fcm] 


Rao [km] 


Al 


1.295 


2.92 


685 


38.00 


13.88 


B1,B2 


1.471 


2.72 


737 


36.93 


13.99 


C1,C2 


1.523 


2.65 


756 


36.71 


14.01 


Dl 


1.649 


2.57 


777 


36.19 


14.03 


D2 


1.649 


2.57 


777 


36.19 


14.02 


El 


1.957 


2.20 


909 


33.40 


13.93 


E2 


1.959 


1.97 


1058 


31.14 


13.07 



the output of the above method, the hydrodynamic quan- 
tities and the conformal factor onto a distribution of SPH 
particles. Finally, this SPH distribution is relaxed to avoid 
spurious inter-particle forces using a braking term 



f = -- 



(37) 



where Virr is the given initial irrotational velocity field of 
the initial mass distribution. Using this method, we produce 
a set of irrotational quasi-equilibrium configurations, both 
varying the EoS and the gravitational mass Moo , where Aiao 
refers to the mass of a single NS in isolation. For Mao, we 
use values between Moo = 1.2Mq which is at the lower limit 
of the observational range and Moo = 1.75Mq which is close 
to the maximal gravitational NS mass of Mmax ~ 1.78Mq 
of the hybrid EoS. The models are summarized in Tab. 
The letter part (A-E) in the model label indicates the mass 
while the number part (1,2) indicates the EoS. All orbits are 
taken slightly outside the ISCO, about 12 ~ 25% in coordi- 
nate separation. This results in smaller lag angles and its os- 
cillations in the initial phase since tidal forces are weaker at 
larger separation. In Tab.|21 we summarize properties which 
result as a consequence of the chosen initial parameters in 
Tab. El 
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3 RESULTS 

In this section, we consider the final evolution of the binary 
NS from a circular orbit configuration to one single object, 
a transient NS or a BH. This evolution takes place on a 
dynamical timescale and is triggered by radiation reaction. 
While in Sec. 18. II the EoS effects on the GW frequency from 
the final inspiral phase just before a merger has been dis- 
cussed, in the latter subsections the result of the dynamical 
evolution of the binary NS from stable circular orbit to the 
merger is presented. 



3.1 EoS impact on the quasi-equilibrium binary 

By considering the compactness {M/R)oa in TabQ we ob- 
serve that the properties of the models become EoS depen- 
dent for the most massive cases Dl, D2 and El, E2, i.e. 
when Moo > 1.5Mq. For lower stellar masses, the NSs do 
not reach the phase transition density pt l.Spo in the 
center and an EoS effect cannot be observed. It has been 
shown that the orbital frequency and the gravitational wave 
frequency of a binar y system depends on the compactness of 
its components (e .g. lLai fc WisemanllQQftlUrvii et al.l200nl : 
iFaber et al.1 1200311 . Therefore we expect that an EoS effect 
may be seen in the GW frequency at the final stage of in- 
spiral in our models. 

Using the method of jUrvu fc Erieuchil l2nont 
lUrvii et al.l l2000l) . we systematically investigated the 
locations of the ISCO, beyond which the binary become 
dynamically unstable, and the GW frequency at the ISCO. 
As shown in Fig. |5Ja), an EoS effect becomes important 
above Moo ~ 1.5Mq. The reason is that at this mass, 
the central density reaches the phase transition density to 
mixed matter in the hybrid case and a dense core of mixed 
matter is formed. For even larger masses, the transition 
density to pure quark matter at p ~ 3.5po is reached and 
a quark core is formed. Due to this, a neutron star with 
the hybrid EoS becomes more compact than one with the 
hadronic EoS as shown in Fig. |Hb), which is reflected in 
the difference in the GW frequency. Near the maximal NS 
mass of the hybrid EoS, the GW frequency of the hybrid 
EoS binary is up to 10% larger than that of the hadronic 
EoS binary. Such a change of tendency for increasing 
gravitational wave frequency with respect to increasing 
mass at Moo ^ 1.5Mq may suggest an existence of a drastic 
change of EoS such as a quark phase transition. 



3.2 Dynamical Evolution of the Neutron Star 
Merger 

We now turn to the merger phase that follows the quasi- 
stationary inspiral when the two NS cross the ISCO and 
become dynamically unstable. To first illustrate the over- 
all dynamics, we plot in Fig. |3 the density distribution for 
model Bl together with the velocity field in the corotating 
frame, whose angular velocity is determined as 



n = 



(38) 



where i runs over all particles and uji denotes the angu- 
lar velocity of the individual particles. The stars are coun- 
terrotating relative to the corotating frame and tidal lag 
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Figure 2. (a)GW frequencies depending on the gravitational 
mass in isolation at the ISCO. (b) Compactness of an isolated 
neutron star with respect to the gravitational mass. 



angles are developing (snapshot 1). At merger, this coun- 
terrotation is first turned into a shear motion at the con- 
tact layer and then into a pattern of vortex rolls due to the 
growing Kelvin-Helmholtz instability (snapshot 2,3). Finally 
we end up with a differentially rotating, non-axisymmetric 
merger remnant (snapshot 4). In Fig. 2] we concentrate on 
the evolution of the merger remnant during the first mil- 
liseconds by plotting the density contour linearly which al- 
lows a better visualisation of the high-density central parts 
(The lowest density contour is at 5 x lO^^g/cm'^.) In this 
figure we can see that during a phase of ~ 2ms after merger, 
the core actually consists of two small subcores which are 
leftovers of the original NS cores. Visible is also that the 
two subcores still carry a small counterrotating motion - 
a direct consequence of the irrotating setup - which disap- 
pears continuously until we end up in an nearly axisym- 
metric, differentially rotating config uration. It is reported 
from finite difference simulations ( e.g. IShibata fc Urvi]l2002l : 
IShibata. Taniguchi fc Urvull2003l) that this twin-core pat- 
tern persists for a longer time than our simulations. This 
may be a consequence of the numerical viscosity in our code. 



3.3 EoS differences in maximal density 

To investigate the different consequences of the EoSs in de- 
tail, we first consider the evolution of the maximum density 
during the merger phase. We split the whole merger event 
in a pre-merger evolution where the binary still consists of 
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Figure 3. Snapshots at selected times of tlie baryonic density distribution and the velocity field of model Bl during pre-merger 
and merger evolution. The density is plotted logarithmically in units of g/cm^, the velocity field is plotted in the corotating frame. 



two tidally stretched objects and in a following post-merger 
evolution where a merger remnant, a NS or BH is forming. 
Since the SPH particle density pa has no direct physical 
meaning, the actual density p{x) has to be calculated as 
a statistical average over the particles densities pa- To find 
the maximal density, we first determine density values on se- 
lected gridpoints and then look for their maximum. Both the 
statistical averaging and the finite grid spacing introduce a 
small noise. In general, the evolution of Pmax, plotted in Fig. 
ISlfor our simulations, shows a slow decrease during the pre- 
merger phase when the two stars become tidally stretched 
and reaches its minimum during the actual merger. This is 
the case when the GW luminosity either approaches its max- 



imum and the two NS are maximally tidally stretched or just 
a bit later when a new merger remnant has already formed 
but a twin-core structure is still present (see sect. l3.2ll . After 
the minimum, we can distinguish three possible evolution- 
ary scenarios of pmax, depending on the collapse behaviour 
of the remnant. If the remnant can be stabilized by pres- 
sure and centrifugal forces, pmax slowly increases and then 
becomes constant. The second scenario is the delayed col- 
lapse where pmax first slowly increases on several dynamical 
timescales. This fragile equilibrium is concluded by a final 
collapse, which is visible as a steep rise of pmax - The last sce- 
nario is the immediate collapse just after the merger. Here, 
pmax increases without delay on a dynamical timescale. 
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Figure 4. Same as Fig|3]but for the post-merger evolution. Tlie density is plotted linearly in units of lO^^g/cm'^. 



Differences in the merger dynamics caused by the EoS, 
appear when pmax crosses the phase transition density pt- 
At this point the hybrid EoS and therefore the dynamical 
evolution of the hybrid models separates from that of the 
hadronic EoS, If the individual stars are more massive than 
Mod — 1,5M0, their central density exceeds the phase tran- 
sition density pt , In this case, EoS differences happen al- 
ready before merger. This makes a very small effect which 
even disappears close to merger due to tidal distortion if Mao 
is just at the threshold of ~ I.SMq as in model D1/D2. On 
the other hand, if Moo is close to the maximal gravitational 
mass of the hybrid EoS as in model E1/E2, the EoS dif- 
ference in Pmax increases to up to ~ 60%, Differences are 
also measurable for global quantities like the angular ve- 
locity of the binary, the compactness and especially for the 
GW frequency (cf, sect, 13. ll . If the initial mass is less than 
Moo — 1.5Mq an EoS difference can only be seen after the 
merger as part of a different evolution scenario of the merger 
remnant. Typically, the remnant collapses immediately or 
after a couple of dynamical timescales in the hybrid cases 



while it settles down to a transient NS in the hadronic cases. 
Among the hadronic models, there is only the very massive 
model El which collapses immediately. In model Dl, the 
maximal density slowly increases but does not really settle 
down - a sign that this object might eventually collapse. 
However, this might well be a numerical effect as angular 
momentum is transported from the core to the outer lay- 
ers and differential rotation is slowly converted into uniform 
rotation. The amount which is transported is on the order 
of 40-50% in the very center during the whole evolution, 
despite of a very good overall angular momentum conser- 
vation, either by numerical viscosity or by gravitational in- 
teraction between the non-axisymmetric core and the outer 
layers. However it is still possible to stabilize a remnant 
with a gravitational mass of nearly SM© with an EoS hav- 
ing a maximal gravitational mass of 2.2Mq for about 12ms 
~ 820A/toi, where Mtot = 2.98M0 is the initial gravitational 
mass for this model. This stabilisatio n effect due to differen- 
tial rotation has been pointed out in JShibata fc Urvi]ll200d: 
iBaumgarte et aLll2000t) . The less massive models C1,B1 and 
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Figure 5. Evolution of the maximal density measured in units of 
the nuclear saturation density po '■= 2-8 X lO^'^g/cm^. The origin 
of the time axis has been shifted to the GW luminosity maximum. 
Solid lines correspond to hadronic models, dashed lines to hybrid 
models. 



Al lead to transient NS which do not collapse on a hydrody- 
namic timescale. On the hybrid-EoS model side, the models 
collapse with different collapse timescales. While models E2 
and D2 collapse immediately, C2 collapses with a delay of 
about ~ 9ms. Model B2 does not collapse within the simu- 
lation time but it also shows a continuously growing maxi- 
mum density. However, these collapse timescales are likely to 
be dominated by the above mentioned angular momentum 
transport. 

3.4 The Gravitational Wave signal 

A second indicator of EoS effects is the GW signal. Since the 
waveform is sensitive to dynamical mass motions we expect 
that all the above maximal density differences are reffected 
in the GW signal mainly in the form of different frequencies. 
In Fig. |S|we plot the waveforms of all models sorted accord- 
ing to the initial mass. Model El /E2, the most massive one is 
the only one which shows a significant pre-merger EoS differ- 
ence in the waveform. This is because the above mentioned 
maximal density difference translates via the compactness 
and the binary angular velocity into a GW frequency differ- 
ence which amounts to ~ 10% at the ISCO (cf. Fig.|5J while 
it disappears during the actual merger phase. The next less 
massive model, D1/D2, does not show any EoS differences 
in the GW frequency at the ISCO and further during the 
pre-merger phase, but we may expect, that for larger binary 
distances and therefore smaller tidal interaction, GW fre- 
quency differences could be seen. The more interesting part 
of this model is the different collapse behaviour and there- 
fore the totally different waveform. While the hybrid model 
D2 only produces a short, high-frequent burst before the 
collapse, the hadronic model emits a long wavetrain which 
decreases slowly in amplitude. Models C1/C2 and B1/B2 
emit both a quasi-periodic (QP) GW-signal which is char- 
acteristic for the rotation and oscillation mode in the merger 
remnant. EoS-differences in the waveform do not appear un- 
til a considerable mixed matter core has formed in the hybrid 
case. As a consequence, the first part of the post-merger GW 
signal, which is the strongest in our simulations, will not be 
affected by any EoS difference. However, when the mixed 
matter core becomes large enough, an accumulating phase 
shift in the waveform becomes clearly visible. In the C1/C2- 
model the shift adds up to more than half a period before 
the collapse of the hybrid remnant happens, in the B1/B2 
model, the shift is only very small, but visible. 

3.5 Gravitational Wave Spectra 

The Fourier spectrum of the GW waveform is plotted in 
FigQ The most massive El/E2-model shows its EoS dif- 
ferences in the pre-merger part, i.e in the frequency do- 
main around IkHz. Here the spectrum produced by the pre- 
merger hybrid waveform is significantly stronger than the 
one of the hadronic waveform, since the amplitude of wave 
strain h scales roughly as {r/M)h ~ {M/R)oo where r is a 
distance from a source. This is consist ent with simulations 
using quasi-equilibrium sequences iFaber et a l. 2003') . There 
is no contribution on the high frequency side from this model 
since both remnants collapse immediately after merger. Note 
that we miss the waveform emitted by the ringing of the re- 
sulting BH. However, the expected frequency of this signal 
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Figure 6. Gravitational waveforms of all models sorted according 
to their initial mass. The origin of the time axis has been shifted to 
the GW luminosity maximum. Solid lines correspond to hadronic 
models, dashed lines to hybrid models. 



lies in tlie ranee of 5 - lOkHz JShibata fc Urvul2002l : Leave j 
llQSSi) we can clearly separate out the BH waveform. 

The other mode[s show their EoS-differences via dif- 
ferent frequency peaks which result from the QP-waveform 
emitted in the post-merger. The difference is most obvious 
in the Dl/D2-mode[. For the hadronic EoS model Dl, a 
distinct Fourier peak is c[early visible, whereas the hybrid 
EoS model D2 [eads only to a much weaker and very broad 
Fourier peak since the remnant is collapsing very soon after 
merging producing a much shorter post-merger GW signal. 
In the C1/C2- model, both the hadronic and the hybrid mod- 
e[s have strong QP peaks where the one of the C2-model is 
sfightly shifted to higher frequencies due to the more com- 
pact remnant. The same can be seen for the model B1/B2. 
However, effects here are smaller and a measurement may be 
difficult. The EoS-effects are this tiny because the spectrum 
is dominated by the first GW burst after merging which is 
insensible to EoS-differences. 

We may compare the obtained QP frequencies . fop to 
value s from a fully relativistic calculation ( Sliibata fc Urvtil 
1^02^ where an irrotational initial velocity field with a poly- 
tropic EoS is considered. Two peaks are found in their spec- 
trum 



/qpi ~ 1-8/qs 
/qp2 ~ 2.8/qs 

for an adiabatic index of F = 2.25 and 



1.8/qs 
3.0/qb 



(39) 
(40) 



(41) 
(42) 



for F = 2. Here, fpE is the GW frequency of the binary at 
the ISCO. dShibata fc Urvii2002^ point out, that fQP2/fQE 
is not very dependent on the initial compactness of the stars 
but the higher peak fQP2/ Jqe depends on F in a way that 
softer EoSs lead to higher values for fQP2/ /qe- Since we 
use realistic EoSs with a non-constant F in our simulations, 
we expec t the ratio fop2/fQE to vary with inc reasing com- 
pactn ess (|Zana,_Caitrella &: McMillan| [l996: Fa ber fc Rasid 
l200lHOechslin et alJl2002^ ■ Indeed our values for fQP2/fQE 
are strongly dependent on the initial mass Moo and com- 
pactness as can be seen in Tab. |3 From the behaviour of 
fQP2/ /qe we can deduce that the EoS softens for more com- 
pact models, i.e. the EoS softens at higher densities. It might 
be very interesting to calculate the fQP2/ fQE-{M / R)oo- 
relation for various other realistic EoSs on the market. Mea- 
suring this relation might give a detailed insight to the NS 
EoS in the nuclear regime above 10^''g/cm^. 

One might expect that /qp2 is overestimated in our 
simulations due the already mentioned numerical viscosity 
effects in the code which leads to a slowly contracting merger 
remnant and therefore to emission of GW with higher fre- 
quencies. To check that this is not the case, we have also 
calculated the spectra of the truncated wavesignals which 
only contain the first few oscillations of the post-merger sig- 
nals. The resulting QP peaks are reduced in strength but 
their positions remains unchanged within the uncertainty of 
the peak width (see Fig. |7J . 
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Table 3. Ratio of /qp2 to Jqe for all models where a transient 
NS forms or a delayed collapse happens. The values for Jqe a^re 
determined from quasi-equilibrium models and are taken at the 
ISCO. 



Model M^[Mq] C={M/R)^ /QP2[kHz] /qp2//qb 



Al 


1.2 


0.128 


1.98 


2.47 


B1/B2 


1.35 


0.143 


2.20 


2.56 


C1/C2 


1.4 


0.148 


2.32 


2.65 


Dl 


1.5 


0.158 


2.70 


2.97 



3.1 




2 4' ' ^ ' ' 
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C 



Figure 8. Relation between the ratio fciPllfQE and the com- 
pactness C = {M/R)ac- The diamonds correspond to the 
hadronic EoS-m odels Al to Dl, the va lues for F = 2 and F = 2.25 
are taken from iShibata fc Urvi^l2002^ 

4 CONCLUSIONS 

In this work, we have considered the impact of the EoS 
to the inspiral and merger dynamics of a binary neutron 
star coalescence. We have performed equilibrium sequence 
studies to investigate the GW frequencies around the ISCO. 
Then, dynamical simulations of the merging event using a 
3-dimensional SPH code have been carried out. We have 
chosen a hadronic EoS based on the relativistic mean field 
approximation and a hybrid EoS which is obtained by 
combining the hadronic EoS with a MIT-bag-model-EoS in 
the high-density regime above pt = l.Spo- The gravitational 
wave frequency at the ISCO depends strongly on the initial 
mass of the binary. For a generic mass of 2 x 1.4A/0, we 
find a GW frequency of around 900Hz for both EoSs. 
EoS differences become important for models with initial 
gravitational masses larger than Mac '~ l.SM©. At this 
mass range, the central density of the individual stars reach 
the phase transition pt included in the hybrid EoS model. 
At masses close to the maximum gravitational mass of the 
hybrid EoS, M = 1.78Mq, the relative difference in faw 
becomes as large as ~ 10%. 

The maximum density evolution during the merger de- 
pends sensitively on the fact whether the phase transition 
density is crossed and at which time. For very massive mod- 
els like models E1/E2 with an initial mass of Moo = 1.75Mq 
this is already fulfilled by the individual companion stars 
in the pre-merger phase. Less massive models with masses 
below Moo ^ 1.5Mq cross this threshold later in the post- 
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merger phase when the matter is contracting to form a tran- 
sient NS or a BH. The hadronic EoS models all form a tran- 
sient NS except the very massive El model which collapses 
immediately. In the Dl model, the maximum density does 
not converge after 12ms but rises slowly and continuously. 
We think that this is an effect of a considerable angular mo- 
mentum transport from the remnant core to the outer lay- 
ers which implies a continuous conversion of the differential 
rotation pattern into uniform rotation. All hybrid EoS mod- 
els collapse on timescales which is highly dependent on the 
model mass. While the massive models E2 and D2 collapse 
immediately, models C2 and B2 collapse after a contraction 
of the remnant on several dynamical timescales and a con- 
siderable decrease of the degree of differential rotation. This 
indicates that the collapse is driven in this case by angu- 
lar momentum transport caused e.g. by numerical viscosity. 
The collapse dynamics of these models will be investigated in 
future work. The emitted gravitational waveforms are very 
dependent on the initial mass of the model. The most mas- 
sive model E1/E2 shows a large difference prior to merger 
both in frequency and amplitude. This is a consequence of 
the different compactnesses for the models El and E2 due 
to the different EoSs. The same effect is not visible in the 
less massive models because the compactness difference due 
to the EoSs vanishes. A second aspect of the GW signal is 
the waveform emitted by the remnant. If the remnant does 
not collapse immediately to a BH, a quasi-periodic wave- 
form with a frequency of 2-3kHz is emitted. The ratio of 
this frequency to the GW frequency at the ISCO fQP2/ fqs 
depends sensitively on the model and therefore on the mass 
or compactness [M/ R)oo- Since fQP2/ fqs is fairly constant 
for an EoS with constant stiffness, we interpret this result 
as a consequence of the varying stiffness of our EoS used. 
Hence, the relation fQP2/ fQE-{M / R)ao is characteristic for 
the behaviour of the EoS's stiffness. A measurement of this 
quantity would provide important additional information of 
the EoS. 

It has been suggested in iHuehed i2002h that measure- 
ments of high frequency GW will require a network of broad- 
band detectors combined with narrowband detectors that 
have a good sensitivity in the high frequency domain. Such 
measurements of high frequency GW may become feasible 
in the future and discover or constrain the existence of phase 
transition in neutron star matter at high densities. 
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